#!/usr/bin/env python

#aim: input stereoisomer smiles
#transform to ---SCI3
#generate 3D conformation
#transform ---SCI3 -> SC

#if unique smile => has to be kekulised (no, that is not anylonger true, I fixed the "unique smiles" converter)

from openeye.oechem import *
from openeye.oeomega import *
import chemistry, mysql
import sys, string




#-----------
if len(sys.argv) <> 4:
	print 'reaction.py  <input_file> <output_file.mol2> <failed smiles>'
	sys.exit(1)

debug = True
ifs = oemolistream(sys.argv[1])
ofs = oemolostream(sys.argv[2])
out_file = open(sys.argv[3], 'w')

#--------------------------------------------
def get_rules():
	cursor = conn.cursor ()
	command = "SELECT smirks, comment, id FROM definition_reactions order by id"
	cursor.execute(command)
	rows = cursor.fetchall ()

	rules = []
	for row in rows:
		rules.append([[row[0]],row[1]])
	
	conn.close()

	return rules
#--------------------------------------------
def reformat_rules(rules):
  #OPENEYE has got a problem with writing [NH3+], ...
  rule_counter = 0
  for rule in rules:
	trans_counter = 0
	for transformation in rule[0]:
		educt, product = transformation.split('>>')
		#primary amines
		product = product.replace('[NH3+]', '[N+]([H])([H])[H]')
		#secundary amines
		product = product.replace('[NH2+]', '[N+]([H])([H])')
		#tertary amines
		product = product.replace('[NH+]', '[N+]([H])')
		product = product.replace('[nH+]', '[n+]([H])')
		product = product.replace('[nH]', '[n]([H])')
		reaction = educt + '>>' + product
		#print reaction
		rules[rule_counter][0][trans_counter] = reaction
		trans_counter = trans_counter + 1
	rule_counter = rule_counter + 1
  return rules
#--------------------------------------------
def fix_bonds (smiles,debug):
#openeye does not write  [NH3+], ... => backward reaction does not work
	old_smiles = smiles

	#fix double bond
    	smiles = smiles.replace('c(=c','c(c')
    	smiles = smiles.replace('c1=c','c1c')
    	smiles = smiles.replace('c=c','cc')
    	smiles = smiles.replace('=c[n','c[n')
    	smiles = smiles.replace('c=n','cn')
    	smiles = smiles.replace('n=c','nc')
    	smiles = smiles.replace('c=[n','c[n')

	if debug and old_smiles <> smiles:
		print 'cleaned smiles: '  + old_smiles + ' -> '+ smiles 
		print
	

	return smiles

#--------------------------------------------------

conn=mysql.connect2server('', 'webuser','purchasable')  			  

mol = OEGraphMol()

iso = True #keep stereoinformation
kek = False #kekulize to get unique smiles (no, this is no longer true)
chiral = True #keep information about chirality (could also be false, because this information is guessed anyway)

reaction = 0
comment = 1

rules = get_rules()
rules = reformat_rules(rules)


result_list = [] #stereoismer smiles generate E/Z -> can result in identical smiles after reaction
#print rules


print 'done'

title_list = [] #this can be placed further down to speed things up but then the user has to ensure that there are no no duplicate molecule names in the list

while OEReadMolecule(ifs, mol):

  #remove spaces in title
  title = mol.GetTitle()
  title = title.replace(' ', '')
  mol.SetTitle(title)
  print title
  #OE3DToInternalStereo(mol)
  OEFindRingAtomsAndBonds(mol)
  OEAssignAromaticFlags(mol)
  OEAssignFormalCharges(mol)
  smiles_list = []
  exists = False #check if a smiles was added to results_list
  smiles = chemistry.CanSmi(mol,iso,kek,chiral)
  old_smi = smiles
  tmp_list = []
  for rule in rules:
    #print rule
    for transformation in rule[0]:
		#make mol
		work_mol = OEGraphMol()
		OEParseSmiles(work_mol, smiles)


		#generate results for all posible matches
		lg = OELibraryGen()
		lg.Init(transformation)
		lg.SetExplicitHydrogens(True)
		lg.SetValenceCorrection(True)
		matches = lg.AddStartingMaterial(work_mol, 0, False)
		if matches > 0 and debug:
			print "number of matches =", matches		
		for product in lg.GetProducts():

			smi = chemistry.CanSmi(product,iso,kek, chiral)
			print "product smiles =", smi
			#print smi, '<----------------------'
			error = chemistry.check_smi(smi, old_smi, transformation, rule[comment], mol.GetTitle(), debug)
			#don't know how to get ride of a leaving group of unknown compestition -> strip, reset error
			OETheFunctionFormerlyKnownAsStripSalts(product)
			if error:
				print 'error!!!!!!!!!!! (or leaving group?)'

			if smi not in tmp_list and not smi in result_list and smi <> smiles:
				tmp_list.append(smi)
				result_list.append(smi)
				if debug:
					print 'added smiles'
			else:
				exists = True
				#tmp_list_kek.append(smi_kek)
			#do not keep org smiles in case of error -> no reaction -> no matching in FlexX

  
  #update smiles_list
  smiles_list = []
  
  if len(tmp_list) == 0 and not exists:
	#did not react
	out_file.write(smiles + '\t' + title + '\n')

  for i in tmp_list:
	smiles_list.append(i)


  
  if debug:
	print 'results'
  	print smiles_list
  
  for moli in smiles_list:

	smi_not_kek = chemistry.clean_smiles(moli)
	stereoisomers = []
	#read smiles
	work_mol = OEMol()
	OEParseSmiles(work_mol, smi_not_kek)
  	OEAssignAromaticFlags(work_mol)
  	OEAssignFormalCharges(work_mol)
	OE3DToInternalStereo(work_mol)

	#first generate correct stereoisomers, only if not specified
	for flip_mol in OEFlipper(work_mol):
		omega_mol = OEMol() #OEGraph does not work
		flip_smiles = chemistry.CanSmi(flip_mol,iso,kek, chiral)
		OEParseSmiles(omega_mol, flip_smiles) #don't now why, but can not run omega on flip mol directly
		#print chemistry.CanSmi(omega_mol,iso,kek, chiral)
		#generate 3D
		omega = OEOmega()
		omega.SetMaxConfs(1)
		#print 'before omega'
		omega(omega_mol)

		#convert SCI3 -> SC

		for atom in omega_mol.GetAtoms():
			atom.SetPartialCharge(atom.GetFormalCharge())
			if atom.GetAtomicNum() == 53:#iodine
				for na in atom.GetAtoms():
					if na.GetAtomicNum() == 6 and not na.IsAromatic():
						omega_mol.DeleteAtom(atom)
						na.SetPartialCharge(-3) #-> easy to recognize for further manipulations

		#generate unique title
		counter = 0
		#print title
		if title in title_list:
			new_title = title + '__0'

			while new_title in title_list:
				counter = counter + 1
				start = new_title.find('__')
				new_title = new_title[:start] + '__' + str(counter)
			title_list.append(new_title)
		else:

			title_list.append(title)
			new_title = title

		#print new_title , '<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'
		
					

		omega_mol.SetTitle(new_title)
		OEWriteMolecule(ofs, omega_mol)
	
out_file.close()
print 'finished reaction'

